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ABSTRACT 

The thermal instabihty with a piecewise power law cooling function is investigated using one- and 
three-dimensional simulations with periodic and shearing-periodic boundary conditions in the presence 
of constant thermal diffusion and kinematic viscosity coefficients. Consistent with earlier findings, the 
flow behavior depends on the average density, (p). When {p) is in the range (1-5) x 10~^''gcm~'^ 
the system is unstable and segregates into cool and warm phases with temperatures of roughly 100 
and 10"* K, respectively. However, in all cases the resulting average pressure {p) is independent of 
(p) and just a little above the minimum value. For a constant heating rate of 0.015 ergg~^ s~^, 
the mean pressure is around 24 x lO^^^dyn (corresponding to p/k-e, ~ 1750 Kcm~'^). Cool patches 
tend to coalesce into bigger ones. In all cases investigated there is no sustained turbulence, which 
is in agreement with earlier results. Simulations in which turbulence is driven by a body force show 
that when rms velocities of between 10 and 30 km/s are obtained, the resulting dissipation rates 
rates are comparable to the thermal energy input rate. The resulting mean pressures are then about 
30 X 10~^^dyn, corresponding to p/k-g « 2170 Kcm^'^. This is comparable to the value expected 
for the Galaxy. Differential rotation tends to make the flow two-dimensional, that is, uniform in the 
streamwise direction, but this does not lead to instability. 
Subject headings: hydrodynamics, instabilities, turbulence, ISM: general 



1. INTRODUCTION 

The importance of thermal instability (TI) has been 
extensively studied in the context of the generation and 
regulation of structures in the atomic interstellar medium 
(the so-called cold and warm atomic phases usually de- 
noted as CNM and WNM), for a review see e.g. Cox 
(2005). The core idea was presented by Field et al. (1969, 
hereafter FGH) in their famous two-phase model: two 
thermally stable phases (cold and cloudy; warm and dif- 
fuse) co-exist in pressure equilibrium regulated by the 
presence of a thermally unstable phase at an intermedi- 
ate temperature. After the observational determination 
of the existence of significant amounts of hot gas in the 
Galaxy, the FGH model was complemented with a third, 
hot, phase by McKee & Ostriker (1977), in which model 
most of interstellar space was occupied by million-kelvin 
gas produced in supernova explosions. Since then, the 
estimates of the filling factor of the hot component have 
reduced to 10%-30% near the Galactic midplane, being 
larger at larger heights. Moreover, most of the hot gas 
seems to be confined in large bubbles created by clus- 
tered supernova activity rather being distributed homo- 
geneously around the Galaxy. In this light, therefore, it 
seems justified to neglect the hot component and return 
to the simpler FGH picture when modeling the colder 
and denser phases of the interstellar medium (ISM). 

There have been a number of numerical investigations 
of the interaction of turbulence and TI. In most papers 
the turbulence is forced by sources other than the TI 
itself: random turbulent forcing at varying scales and 
Mach numbers (e.g. Gazol et al. 2005), localized injec- 
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tions of energy mimicking stellar winds (e.g. Vazquez- 
Semadeni et al. 2000), the magnetorotational instability 
(e.g. Piontek & Ostriker 2005), and systematic large- 
scale motions such as propagating shock fronts (e.g. 
Koyama & Inutsuka 2002) and converging flows (Au- 
dit & Hennebelle 2005) have been considered. One of 
the major findings from these models is that, because 
of the turbulence present in the system, large pressure 
deviations are generated and significant amounts of gas 
can exist in the thermally unstable regime. These re- 
sults suggest that the FGH picture of the ISM exhibiting 
"discrete" temperatures and densities and a unique equi- 
librium pressure should be modified in the direction of a 
"continuum" of states with an overall pressure balance 
but with large deviations from it. 

In recent years the possibility of driving turbulence by 
the TI itself has received some revived attention. Con- 
trary to Kritsuk & Norman (2002a), who found turbu- 
lence to die out as a power law, Koyama & Inutsuka 
(2006) find the turbulence to be sustained — at least for 
times up to 0.1 Gyr. The possibility of TI- induced tur- 
bulence is potentially similar to the Jeans instability in 
a self-gravitating medium that is able to maintain a sta- 
tistically steady state in which the instability drives the 
turbulence and turbulent heating prevents the disk from 
cooling into a static equilibrium. Simulations by Gam- 
mie (2001) have shown that such a state of self-sustained 
"gravito-turbulence" is indeed possible. Wada et al. 
(2000) found a similar result for the case of the com- 
bined action of gravitational and thermal instabilities. 
The possibility of driving turbulence by means of insta- 
bilities is indeed quite common in astrophysics. Espe- 
cially popular is the magneto-rotational instability that 
is known to drive turbulence in disks (Hawley et al. 1995, 
Brandenburg et al. 1995), but there is also the Rayleigh- 
Benard instability, which leads to turbulent convection 
(e.g., Kerr 1996). 



In this paper we focus on the interaction between 
turbulence and the nonhnear stages of the TI, starting 
from one-dimensional calculations and extending them 
to three dimensions. Following an approach similar to 
those of Koyama & Inutsuka (2004) and Piontck & Os- 
triker (2005), we include thermal conduction, which sta- 
bilizes the gas at wavelengths smaller than the criti- 
cal wavenumber of the condensation mode (Field 1965). 
This wavelength is usually referred to as the the Field 
length; it allows the structures generated by the TI to be 
resolved by the chosen numerical grid. Other approaches 
have also been used: In the model of Sanchez-Salcedo 
et al. (2002), a nonuniform grid was used to resolve all 
the scales down to the cooling length, but nevertheless 
the required amount of gridpoints restricted the calcu- 
lations to one dimension. In some models (e.g., Gazol 
et al. 2005), no bulk viscosity or thermal conduction is 
used, but they are replaced by local resolution-dependent 
artificial viscosities damping Nyquist-scale unresolvable 
structures. 

2. MODEL 

2.1. Governing equations 

We consider the governing equations for a compressible 
perfect gas, 

' =-V-u, (1) 
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where u is the velocity, p is the density, s is the specific 
entropy, with S^ = ■^{uij + Uj^i) — ^Sij'V ■ u being the 
traceless rate of strain tensor, v is the kinematic viscos- 
ity, X is the thermal diffusivity, and C is the net cooling 
or heating, that is, the difference between cooling and 
heating functions, with 

£ = pA - F, (4) 

where F = const is assumed for the heating function. 
Here we consider the photoelectric heating by interstellar 
grains caused by the stellar UV radiation field, for which 
Wolfire et al. (1995) give a value of 0.015ergg^^ s~^ at 
n — 1 cm~^. 

Following common practice, we adopt a perfect gas for 
which p and s are related to pressure p and temperature 
T by the relations 
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where TZ — 8.314 x 10^ cm^ s^^ K^^ is the universal gas 
constant, p is the mean molecular weight (here we as- 
sume /i = 0.62 in all cases and neglect the effects of 
partial ionization), and TZ/p = Cp — Cy, with Cp and Cy 
being the specific heats at constant pressure and volume, 
respectively; 7 = Cp/cy = 5/3 is their assumed ratio. 
The adiabatic sound speed Cg and the temperature are 
related to the other quantities via c^ = jTZT/p. The spe- 
cific entropy is defined up to a constant sq, whose value 
is unimportant for the dynamics. 

We adopt a parameterization of the cooling function 
approximately equal to that given by Sanchez-Salcedo et 



TABLE 1 
Coefficients for the cooling curve given by eq. (6). 
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al. (2002), which has been obtained by fitting a piecewise 
power law function of the form 

A{T) ^ Cu+iT^-^+' for T, < T < T,+i, (6) 

to the equilibrium pressure curve of the standard model 
of Wolfire et al. (1995) for the ISM in the solar neighbor- 
hood. When thermal equilibrium with the chosen back- 
ground heating function, F is assumed, and a continuity 
requirement for A, 



C,- 



-l.i 



-a 



iS+l 



T, 



(ft,.+i-/3._i,i) 



(7) 



is taken into account, we arrive at the values of the co- 
efficients listed in Table 1. The coefficients Ci^i+i given 
by Sanchez-Salcedo et al. (2002) deviate from this con- 
dition by 4%-8%. It turned out that with their original 
coefficients the flow amplitude showed spurious oscilla- 
tions in time which disappeared when we use the revised 
coefficients. 

It is convenient to measure time in gigayears, speed 
in kilometers per second, and density in units of 
10"^'* gcm~'^. Pressure is therefore measured in units 
of lO^^^dyn. Our unit of length is therefore l(km/s) x 
1 Gyr = 1.02 kpc; in the following, we denote the unit 
of length for simplicity as 1 kpc, keeping in mind that it 
should really be 1.02 kpc. Viscosity and thermal diffu- 
sivity are measured in units of Gyr km s~^. 

We use periodic boundary conditions in all three di- 
rections for a computational domain of size (200 pc)'^, 
which is the typical domain size employed in simula- 
tions of supernova-driven turbulence in the interstellar 
medium. However, smaller domains would be more suit- 
able to resolve the smaller scales, as has been done by 
Kritsuk & Norman (2002a), for example. We use the 
Pencil Code,^ which is a non-conservative, high-order, 
finite-difference code (sixth order in space and third or- 
der in time) for solving the compressible hydrodynamic 
equations. Because of the non-conservative nature of 
the code, diagnostics giving the total mass and total 
energy (accounting for heating/cooling terms) are mon- 
itored and simulations are only deemed useful if these 
quantities are in fact conserved to reasonable precision. 
The mesh spacings in the three directions are assumed 
to be the same, that is, 6x = 6y = Sz. 

We emphasize that no shock or hyperviscosity has 
been used in the present simulation. Therefore, the only 
means of stabilizing the code is through regular viscosity 
v and thermal diffusivity x- In order to damp unresolved 
ripples at the mesh scale Sx in a trail of structures mov- 
ing at speed U, the minimum viscosity and minimum 
diffusion must be on the order of 0.01 USx (see Branden- 
burg & Dobler 2002). In all our simulations the veloc- 
ities are subsonic, so the fastest pattern speed is given 

* see http://www.nordita.dk/softwarc/pcncil-code. 



by the sound speed. In the foUowing we quote the mesh 
Reynolds number based on the mean (volume averaged) 
sound speed, Cg, and the mesh size Sx, 
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CsSx/l^, 
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The minimum viscosity quoted above corresponds to a 
largest permissible value of Rcniosh of about 100. How- 
ever, in the presence of strong converging flows and 
shocks the largest permissible value may be of order 
unity. 

Since we want to use minimal values for v and x in both 
the warm and cold components we keep v and x constant 
rather than, for example, the dynamical viscosity or the 
quantity JC = pjx (see, e.g., Piontek & Ostriker 2004). 
In the latter case, x would vary by 2 orders of magnitude 
between warm and cold phases. If the mesh were suffi- 
ciently fine, one could allow for a physically motivated 
dependence of x on T, but this is neglected here. 

In the calculations, we have adopted two different val- 
ues of 1^ and X (5 X 10^'^ and 5 x 10^'* Gyrkm^ s^^), 
keeping their ratio, the Prandtl number Pr ~ v/Xi fixed 
to unity. The corresponding Field lengths, calculated 
from equation (12) using the initial cooling timescale 
of approximately IMyr, are 24 and 7.7 pc, respectively. 
Compared with the average value of the thermal dif- 
fusion in the neutral ISM, roughly 6 x 10^°cm^s^^ w 
2 X 10~^ Gyrkm s~^, corresponding to a Field length of 
about 0.5 pc, the adopted values are larger by 2-4 or- 
ders of magnitude. The cooling length /cool ~ TcooiWrnis 
is close to the physical Field length, being roughly 0.4 pc. 
Our chosen values of x due to the preference of a large do- 
main size, are therefore too large to resolve the fine struc- 
ture in the accretion fronts that result from the cooling 
process. This is a similar setup to the one investigated by 
Piontek & Ostriker (2004, 2005); models achieving Field 
lengths smaller than the cooling length include for ex- 
ample, Sanchez-Salcedo et al. (2001), Kritsuk & Norman 
(2004), and Koyama & Inutsuka (2004). 

2.2. Stability properties 

The first thorough stability analysis was done by Field 
(1965), who also included the stabilizing effect of thermal 
diffusion. Assuming the solutions to be proportional to 
exp(nt + ifc • a;), the dispersion relation can be written in 
the form 



n{n + n^){n + l3np + n^)+ujl 
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-0, 

where we have also included the effect of kinematic vis- 
cosity. Here, Wac = Cgfc is the acoustic frequency and 
/3 — dlnA/dlnT is the local logarithmic slope of the 
cooling function. We have restricted ourselves to cases 
where F is constant and A depends only on T. The cool- 
ing time is characterized by the quantity 



np = poCp/{cyT), 



(10) 



which is to be evaluated for the equilibrium solution. 



Here, 



'Cp - 



{dC/dp)T = A. Note that Up is just 



the inverse cooling time defined by Piontek & Ostriker 
(2004). The subscript p follows from a similar notation 
used by Field (1965) who defined instead a wavenum- 
ber kp — rip/cs, which is also referred to as the cooling 
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Fig. 1. — Dispersion relation n{k) in the unstable regime with 
/3 = 0.56, obtained by solving eq. (9) for a representative range 
of values of np/(csA;p). The n(fc) curves are normalized in terms 
of fcp and Cgfep. In the range np/{cskY) = 0.2-2 (solid lines) the 
maxima of n/(cBkp) are monotonically increasing. The curves for 
np/{csk-p) = 5 (dashed line) and np/(csfcF) = 10 (dotted line) 
deviate from this trend. The diagonal dash-dotted line indicates 
the approximation valid for small wavcnumbcrs [eq. (13)]. 

wavenumber. Viscous and diffusive effects are character- 
ized by the corresponding rates. 
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Thermal instability is only possible for /3 < 1. This con- 
dition corresponds to the isobaric instability criterion of 
Field (1965). The isochoric and isentropic criteria, /? < 
and /3 < —1/(7 — 1) = —3/2, respectively, are less strict 
in that the isobaric criterion for instability is then auto- 
matically satisfied. 

When thermal diffusivity is included, the gas can be 
stabilized (even though /3 < 1) provided the largest pos- 
sible wavenumber in the system (which we denote as ki ) 
is larger than the Field wavenumber, k-p, defined as 



(1 - /?)np/(7x) (for /3 < 1). 



(12) 



The instability has therefore the character of an ordinary 
long- wave instability requiring ki < kp. The correspond- 
ing dispersion relation is shown in Figure 1 for various 
values of Up using ly = x (top) and ly ^ (bottom). 
The value of x is given in terms of the ratio np/{cskp), 
for which three values have been chosen to illustrate this 
dependence. As expected, the presence of kinematic vis- 
cosity has a stabilizing effect. Setting v — leads to 
somewhat larger growth rates, especially when np/{cskp) 
is large and n/{cskp) is small. For np/{cskp) = 2, for 
example, the normalized growth rate for v = x ^^ the 
largest among the three cases shown in Figure 1, but it 
hardly increases when i^ — > 0, in which case the growth 
rate is actually the smallest among the three cases. 

In the limit fc ^ fcp , which is relevant when diffusive ef- 
fects are negligible, the unstable branch of the dispersion 
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Fig. 2. — Net cooling vs. temperature for three values of p, given 
in units of \p] = 10"^^'* dyn. 

relation reduces to 



== Cisofc V/3 



1, 



(13) 



where we have introduced the isothermal sound speed, 
Ciso = Cs/^/j- Note that in the thermally stable case 
with /? 3> 1 we obtain the usual dispersion relation for 
isothermal sound waves, to = Cisok, where a; = in is the 
frequency. For (3 = 1/2 this approximation yields n = 
Cisok, as stated by Field (1965) in his equation (36). For 
P = 0.56 this approximation is shown in Figure 1 as a 
straight dash-dotted line. 

2.3. Saturation properties 

In the absence of thermal diffusion, thermal equilib- 
rium is given by the condition C = 0. Pressure equi- 
librium between the cold and warm phases requires that 
equilibrium is achieved under the constraint of constant 
pressure. Such an equilibrium would however only be sta- 
ble if a temperature increase would lead to correspond- 
ingly more cooling, that is, if 



dC 
df 



— >0, (stability). 



(14) 



where the subscript p indicates that the pressure is held 
constant. In Figure 2 we plot £ as a function of T for 
constant p; three values are considered: p = 25, 35, and 
50, all in units of \p] = 10~^^ dyn. This figure shows that 
there can be two stable states at about 10^ and 10^ K. 
We denote these values by Tq and Tyj for the cold and 
warm phases. At T « 10'^ K there is an unstable equilib- 
rium, whose temperature is denoted Ty. The densities of 
the three equilibria, obtained by solving C{T;p) = for 
T numerically for given p and then expressing the result 
in terms oi p = p{T,p), are plotted in Figure 3. 

When setting up a simulation the density is particu- 
larly useful, because its mean value in a certain volume 
is proportional to the mass, which is constant for closed 
and periodic boundary conditions, such as those consid- 
ered here. Thus, one can ask the question what is the 
resulting mean pressure as a function of the mean den- 
sity. Of course, as long as the gas is thermally stable, the 
density will be uniform and hence its mean value is al- 
ways equal to the actual value at any point, so it is given 
by combining the equation of state with the condition of 
thermal equilibrium. As is evident from Figure 3, when 
the density is in the range (1-5) x lO^^^gcm^, there is no 
stable solution. This means that the gas will fragment 



Fig. 3. — The two stable solution branches, pc ^ind pyj (solid 
lines), and the unstable solution branch, pu (dotted line), as a func- 
tion of p. On the top the pressure is normalized by the Boltzmann 
constant, p/k-Q. 

into cold patches of temperature T = Tq with density 
PC I a-iid the rest of the ambient gas warms up to the sta- 
ble solution branch T = Tw with density pw- As a direct 
result of mass conservation in our periodic domain, the 
filling factor of the cold component can be expressed in 
terms of the mean density, (p) , which is known from the 
initial condition. Using the definition of the filling factor. 



/pc + (l-/)pw = (p), 
the value of / is given by 
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PC - Pw 

A similar analysis can also be adopted for calculating 
(T) . This allows us to calculate the correlation coefficient 
e in the relation 



{pT)=e{p){T) 



where 



T, 
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0.013 
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The correlation coefficient is small because {pT) de- 
creases slightly and (p) (T) increases strongly as the sys- 
tem segregates, as demonstrated below (in connection 
with Figure 4). The expression {p){T) is almost entirely 
determined by the product of the volume-weighted den- 
sity (or relative mass) in the cold phase, /pc, and the vol- 
ume weighted temperature in the warm phase, {I — f)Tyf, 
so both factors are large compared with their respective 
average values. 

The segregation phenomenon has already been stud- 
ied in a one-dimensional model (Sanchcz-Salcedo et al. 
2002). Here, except for an additional perturbation, the 
initial condition is assumed uniform, p = po = (p), and 
the value of po is varied between different simulations. 
In all the runs presented below the initial perturbation 
is gaussian noise with an rms fluctuation amplitude of 



When Po (in units of 10 gem ^) is be 



10^26 g(,jj^-3 

tween 0.96 and 5.1, the gas is thermally unstable and seg- 
regates into cold and warm components. As time goes on, 
some of the cold spots may move because of slight pres- 
sure imbalance until they coalesce into bigger fragments. 
This coalescence was also found by Sanchez-Salcedo et 
al. (2002) and Koyama & Inutsuka (2004). 

In Figure 4 we plot the evolution of InT in a space- 
time diagram (top) and that of the mean pressure in 
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Fig. 4. — Evolution of In T in a space-time diagram (top) and of 
the mean pressure (bottom) in a one-dimensional simulation with 
1024 mesh points and u = x = 5 X 10~* Gyr km s"^. During 
early times, the rms velocity grows exponentially at a rate of about 
220 Gyr" 1. 

a one-dimensional simulation. Here, v — x — ^ ^ 
10^'* Gyr km^ s^^ which, together with the initial val- 
ues of Cs — 7.5kni/s and Up = 980Gyr~ , yields fcp = 
720 kpc^^ = 23fci, and hence np/(csfcF) ~ 0.2. 

During early times the rms velocity grows exponen- 
tially at a rate of about 210Gyr~^, which is consistent 
with the peak value of n/{csk-p) ~ 0.04 for our set of 
parameters. Note that the mean pressure settles around 
24 X 10~^^dyn once the instability has saturated. At 
that time, smaller structures may still coalesce into larger 
ones, but the total filling factor remains approximately 
constant. During the evolution away from the unstable 
homogeneous state, the mean pressure (proportional to 
{pT)) decreases by about a factor of 2, but the product 
{p) (T) increases by almost a factor of 4. When po is be- 
tween 5.2 and 11 (in units of 10~^''gcm^^) the gas is 
marginally stable (/3 = 1; see Table 1), so in that range 
there will be no segregation into different phases. 

When the mean density is outside the range between 
0.96 and 5.1 (in units of 10~^''gcm~^), the gas is ther- 
mally stable and remains uniform. The dependence of 
the pressure on the density can be obtained in parametric 
form by calculating, using temperature as a parameter, 
p{T) andp(T), that is, 

,(T).^, P(T)^5^^, (19) 

and plotting the two against each other (see Figure 5, 
dotted line). The numerically obtained values for the 
mean pressure {p) , for different mean densities (p) , agree 
with those obtained under the assumption of homogene- 
ity. 

In the unstable regime the pressure is, surprisingly, 
independent of (p), and always around {p) « 24.2 x 
10"'^'' dyn. (Figure 5 shows slight variations about this 
value; this is probably a consequence of the fact that the 
coefficients C^^i+i are only implemented up to three to 
four significant figures, so the cooling curve is still not 
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Fig. 5. — Mean pressure vs. mean density in a one-dimensional 
simulation (solid line), compared with the values obtained for a 
homogeneous system (dotted line). The dotted line (which agrees 
with the solid line in the stable regime and hence cannot be seen 
there) was obtained by plotting p(T) vs. p{T) using T as a pa- 
rameter in eq. (19). On the right axis, the pressure is normal- 
ized by the Boltzmann constant, p/ks, and at the top the num- 
ber density is given. The simulation has 128 mesh points and 
u = x = 5 X 10~^ Gyr km^ s~^. 
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Fig. 6. — Filling factor / as a function of the mean density, as 
predicted by eq. (16). 

perfectly continuous.) Figure 3 shows that for (p) ^ 
24.2 X 10"^'* dyn the warm and cool phases have pw ~ 
0.19 and pc ~ 14.3, respectively. This allows us to deter- 
mine the filling factor as a function of (p); see Figure 6. 
In most of the runs considered below we expect (p) =2, 
so / « 13%. In practice we estimate the filling factor 
as the fraction of mesh points for which T < Tu, where 
Tu ~ 420 K (corresponding to pu = 4.3 x 10^^"* g cm"'^ in 
Figure 3) for (p) = 24.2 x 10"" dyn. The fiUing factors 
determined in this way are quoted for the simulations 
presented below. 

There is a tendency for cool patches to travel and to 
coalesce into bigger ones (see, e.g., Sanchez-Salcedo et al. 
2002). This property is reminiscent of earlier work in the 
context of the thermal instability. Elphick et al. (1991, 
1992) found traveling front solutions and also the merg- 
ing of smaller patches into bigger ones, which they asso- 
ciate loosely with an inverse cascade behavior. However, 
in their work they only discuss the energy equation and 
not dynamical processes. In the case they discuss the 
kink and antikink fronts always travel toward or away 
from each other, thus resulting in the annihilation and 
creation of denser clouds. This is not seen in the present 
work. Also, they discuss much smaller objects of size 
~ 0.02 pc which have considerably shorter sound cross- 
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Fig. 7. — The rms velocity and kinetic energy for a three- 
dimensional run with ;/ = x = 5 X 10^"^ Gyr km s~^ and 256"^ 
mesh points. The dashed line represents a solution with !^ = X = 

5 X 10"'* Gyr km s~^, which was restarted at t = 0.1 Gyr from a 
run with 10 times higher viscosity. Here (p) K, 1.7 X 10"^'* gcm~^, 
(p> ^ 24.3 X 10"" dyn, and (T) sd 8200 K. 

ing times. Furthermore, early on in their evolution our 
clouds tend to accelerate toward each other, as can be 
seen from the curved trajectories. 

3. THREE-DIMENSIONAL SIMULATIONS 

3.1. Fully 'periodic boundary conditions 

In this section we discuss the results of three- 
dimensional simulations. The basic properties of the one- 
dimensional simulations, presented in §2.3, carry over to 
the three-dimensional regime. As expected, the growth 
rates are the same as those found in the one-dimensional 
case. The resulting mean pressure (p) and hence the 
filling factor, as given by eq. (16), are also quite simi- 
lar to those of the one-dimensional case. Nevertheless, 
even though significant amounts of turbulent heating are 
being produced at the most violent phase of the instabil- 
ity, there is in our cases always a subsequent relaxation 
phase in which the flow speed tends to vanish on a long 
time scale (see Figure 7). This agrees with earlier find- 
ings of Kritsuk & Norman (2002a) . The energy decay is 
consistent with a t~^'^ law, just like in ordinary turbu- 
lence (e.g. Mac Low et al. 1998, Haugen & Brandenburg 
2004) . This is also consistent with the results of Kritsuk 

6 Norman (2GG2a), who reported decay exponents in the 
range 1-2 for box sizes between 5 and 500 pc using also 
a more detailed cooling curve in tabular form. On the 
other hand, Koyama & Inutsuka (2006) find that turbu- 
lence remains self-sustained for times up to 0.1 Gyr. 

We emphasize again that we have used constant kine- 
matic viscosity and constant thermal diffusivity in our 
simulations. For the runs shown in Figs 7-13 we have 
used ly — X — 5 X 10^^ Gyr km^ s~^ until t = 0.1 Gyr 
(corresponding to Rcmosh = 2). This corresponds to 
fcp — 230 kpc^ = 23fci, and hence np/{cskp) w 0.6, so 
the initial growth rate is 160 Gyr^ . This is again con- 
sistent with Figure 4 yielding a peak value of n/{csk-p) ss 
0.09 for our set of parameters. 

However, after having reached the peak velocity, the 
flow has becomes sufficiently quiescent so that it is pos- 
sible to decrease the viscosity by a factor of about 10, 
corresponding to Remesh = 20. Figure 8 shows images 
of In T on the periphery of the simulation domain at a 
few selected times after having lowered the viscosity and 
thermal viscosity. Animations of temperature and den- 



sity^ show that late in the simulation, cold patches of gas 
are still moving about, but this is presumably just a re- 
sponse to small-amplitude, small wavenumber variations 
in overall pressure requiring a much longer time scale to 
equilibrate. 

During the course of the simulation the value of fcp 
(based on the averaged value of rip) increases between the 
initial value before saturation of the instability {kpSx « 
0.2) and the saturated state {k-p5x « 1 with the higher 
viscosity and (k-pdx « 3 with the lower viscosity). At the 
end of the simulation the gas is sharply segregated into 
warm and cool phases in almost perfect pressure equilib- 
rium. This can be seen clearly from probability density 
functions of the various quantities that are discussed in 
§3.3. 

3.2. Shearing-periodic boundary conditions 

The shearing sheet approximation simulates the local 
conditions in a disk with strong radial differential rota- 
tion in the limit of large radii. Curvature can thus be 
neglected and the shear can be assumed linear in ra- 
dius, so that we only have an underlying linear shear 
flow C/q = (0, S'a;,0) , where S is constant. The Cori- 
olis force, 211 x u. is added, where ft ~ (0, 0, ft) is the 
angular velocity vector. It is assumed that S scales with 
the angular velocity, so here we take S = —fl which is 
appropriate for galactic disks with a constant linear ve- 
locity law. The combined effects of shear and Coriolis 
force can be subsumed into a single vector (Brandenburg 
et al. 1995), 



C2flUy 
-{2n + s)u, 



(20) 



which is then added on the right hand side of eq. (2). 
After this modification the velocity u describes the de- 
viation from the shear flow and does thus not include 
the basic shear. The basic shear flow appears still ex- 
plicitly as an additional advection operator of the form 
C/o • V = Sxdy. 

In the following we consider the case fl = 100 Gyr^^, 
but we have also considered the case ft = 25Gyr~^ (ap- 
propriate for our Galaxy). The difference between the 
two simulations is small. The main thing that happens 
in all these simulations is a tendency for the flow to be- 
come sheared out, so any variations in the streamwise 
direction become sheared out and the flow becomes es- 
sentially two-dimensional; see Figure 9. However, shear 
does not seem to lead to instability, even though the kine- 
matic growth rate of the thermal instability is apparently 
somewhat increased (190 Gyr^^ instead of 160 Gyr^^). 
This absence of sustained turbulence is somewhat disap- 
pointing, because one might have hoped that the thermal 
instability would have led to condensation in the stream- 
wise direction and thus to new structures that could then 
be sheared out again. This seems to be prevented by 
the general tendency of coalescence, preventing breakup 
into new structures in the streamwise direction. How- 
ever, it may still be interesting to reconsider this issue 
in the future at signiflcantly higher resolution and larger 
Reynolds number. 

^ see http://www.nordita.dk/ brandcnb/niovics/thcrniaLinst. 
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np/{cBkp) = 1.5. The growth rate is about 190 Gyr" , which is somewhat larger than for the corresponding non-shearing run. Note that 
the initially produced structures are quickly sheared out. 

3.3. Forced simulations 

Given that the TI did not produce sustained turbulent 
flows, we consider now cases which turbulence is driven 
by an additional body force in the momentum equation. 
We consider here a forcing function consisting of plane 
waves whose wavevector is chosen randomly at each time 
step and has length between 2.5 and 3.5 times the small- 
est wavenumber in the box, fci = 27r/(0.2kpc). This 
forcing function is therefore (5-correlated in time and ap- 
proximately monochromatic in space (see also Sanchez- 
Salcedo et al. [2002] for simulations in one dimension). 

It turns out that when the flow is driven sufficiently 
strongly to produce rms velocities of around 10-30 km/s, 
the turbulent energy that is dissipated into heat is only 
about comparable to the energy needed to balance the 
losses from cooling (see Figure 10). The mean pressure is 
increased slightly to about 30 x 10"'^"' dyn, corresponding 
to p/fce ~ 2170 Kcm~^. In both cases the spectra of 
u (kinetic energy) and p (density) are similar, except 
that the unforced run shows more relative power in the 
density spectra at large scales; see Figure 11. Over a 
small range of wavenumbers the local slope of the kinetic 
energy spectra is around —5/3. By comparison, Kritsuk 
& Norman (2004) found shallower spectra with spectral 
slope close to —1 in their decaying simulations with TI, 
but this could be a feature of the numerical dissipation 
used in their code. The dissipation wavenumber, fed = 
(uji-nis/i^)^''^, where ojrms is the rms of the vorticity, u) — 
V X M is shown for the unforced run. For the forced run 
the adopted viscosity is critically low, as evidenced by 
the small rise in the kinetic energy at large wavenumbers. 
In fact, the dissipation scale for this run is just outside 
the plot range. It is perhaps because of the presence of 
cooling, which contributes to energy removal, that this 
run has still been successful. 

As in the unforced case, the gas is segregated into warm 




Fig. 10. — Comparison of the rms velocity (top), dissipation rate 
EK (middle), and mean pressure (bottom), for a forced simulation 
(solid lines) and a nonforced simulation (dotted lines; same run as 



in Figure 8). Both cases have u ■- 
256^ mesh points. 
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and cool phases, but now they are only in approximate 
pressure equilibrium; in Figure 12, we show probability 
density functions (PDFs) of Inp, InT, and Inp. The tur- 
bulent increase of the mean pressure has only a small 
effect on the preferred temperatures in the warm and 
cold phases, whereas the density peaks are shifted to- 
ward higher densities, as expected if the system were 
still following the equilibrium pressure-density relation. 
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The turbulence forced at relatively small scales has the 
most drastic effect on the cold cloudy component, the dis- 
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Fig. 13. — Two-dimensional probability density functions of Inp 
and Inp for forced (top) and unforced (bottom) simulations. The 
solid line indicates the thermodynamic equilibrium solution. Dark 
shades indicate large values of the probability density. 

tribution of which has become significantly wider while 
the high density wing was developing. The maximum 
density in the forced case is roughly an order of magni- 
tude larger than in the pure TI case. A similar wing is 
observed at low temperatures, reaching values down to 
the cooling cut-off of 10 K in the highest density regions. 
While in the pure TI case the pressure in the saturated 
state shows a very narrow distribution around the mean, 
in the forced cases the distribution is broad with extrema 
that vary by almost 1 order of magnitude. In addition 
to this broadening of the pressure distribution already 
pointed out in several previous studies (e.g. Gazol et al. 
2005), the amount of gas in the "forbidden" (thermally 
unstable) regime has been observed to increase; in our 
calculations this is seen as a systematic increase of the 
level of the PDFs in between the two preferred states, 
while the peaks themselves become less pronounced. In 
the forced case, about 6% of the gas is found in the unsta- 
ble range where p is between 1 and 5 times 10"^"' g cm~^. 
In the pure TI case, on the other hand, only 2% is in this 
range. In addition there is a significant fraction of cold 
high-density overpressured gas that is in the thermally 
stable regime. Nevertheless, even in the highly turbulent 
regime the signatures of pure TI are still clearly visible 
in the density and temperature PDFs (better for warm, 
worse for cold); the pressure PDF develops broad wings, 
as is familiar from supernova-driven turbulence simula- 
tions (e.g., Korpi et al. 1999; Mac Low et al. 2005), and 
some earlier TI simulations with forced turbulence (e.g. 
Gazol et al. 2005). Still, the mean pressure determines 
the preferred densities and temperatures in the warm 
and cold phases as though the system followed the equi- 
librium pressure-density relation. 



It is customary to discuss scatter plots of pressure ver- 
sus density (Sanchez-Salcedo et al. 2002, Piontek & Os- 
triker 2004, 2005), which allow one to discuss the degree 
to which the gas is locally in equilibrium. In Figure 5 we 
showed that the mean pressure, that is, averaged over 
the entire box, is « 24 x 10~^^dyn when the mean den- 
sity is in the unstable range, {p) = (1^5) x lO^^'^gcm"^. 
It turns out that this result also holds locally, as can be 
seen from a scatter plot of pressure versus density and, 
more conveniently, from a two-dimensional PDF (Fig- 
ure 13) showing the logarithm of the probability density 
as a function of both Inp and \np for both forced and 
unforced runs. In the unforced case the local pressure 
is concentrated sharply around 24 x 10~^'*dyn over a 
broad range of local densities, p ~ 0.2-20 x lO^^'^gcm"^. 
In the forced case, the distribution is broadened around 
the average pressure for p = 0.2-20 x 10~^''gcm~^, but 
there are also dense spots with p > 20 x 10~^^gcm~'^ 
where the gas follows the equilibrium distribution quite 
sharply, confirming earlier findings of Sanchez-Salcedo et 
al. (2002) and Piontek & Ostriker (2004, 2005). 

4. CONCLUSIONS 

Our results confirm the basic findings of Kritsuk & 
Norman (2002a) in that the TI does not lead to self- 
sustained turbulence. In the cases considered in this 
paper the instability just leads to segregation into two 
different phases, and produces only small velocities in 
response to the remaining pressure fluctuations. While 
the growth of the instability occurs over relatively short 
timescales of a few tens of millions of years, the kinetic 
energy of these motions decays exponentially with a slope 
consistent with —1.2 leading to insignificant rms veloci- 
ties after a few hundred millions of years. Thus, in agree- 
ment with Kritsuk & Norman (2002a), the TI alone does 
not lead to self-sustained turbulence. This is somewhat 
different from the purely two-dimensional TI cases inves- 
tigated by Piontek & Ostriker (2004), who report weak 
(« 0.5 kms""'^) non-decaying turbulence over timescales 
of roughly 0.5 Gyr. This behavior seems to carry over 
into three dimensions (Piontek & Ostriker 2005, see their 
Fig. 11). 

Similar results have also been found recently by 
Koyama & Inutsuka (2006) , who also include an explicit 
dynamical viscosity. For times up to 0.1 Gyr their results 
are nevertheless qualitatively similar to ours in that they 
also report rms velocities in the range 0. 1-0.4 km/s, and 
their flow topology is similar to ours at early times. They 
also study smaller box sizes, but their highest turbulence 
levels occur for their largest box size of L = 144 pc, which 
is similar to ours. In both cases the Field length is about 
1/100 of the box size. However, if there is really a differ- 
ence in sustaining turbulence over long times, then this 
might be due to a different formulation of thermal con- 
duction which varies here with density, but is constant in 
the simulations of Piontek & Ostriker (2004, 2005) and 
Koyama & Inutsuka (2006). In the latter case a constant 
dynamical viscosity is included, while in our case a con- 
stant kinematic viscosity is used. Another difference is 
the discontinuous nature of the transitions in the previ- 
ously used cooling function. (The latter was observed to 
lead to spurious oscillatory motions in some of our pre- 
liminary investigations that arc not reported here.) How- 
ever, if the turbulence is real, then this could perhaps be 



understood as an analogy to the Tl-driven turbulence 
found by Kritsuk & Norman (2002b) in the presence of a 
time-dependent heating rate. The idea would be that a 
variable heating rate could perhaps be simulated by in- 
troducing nonlinear feedbacks in some of the coefficients. 

Another possibility for driving turbulence has been dis- 
cussed by Murray et al. (1993). They find that a sys- 
tem segregated into two phases by the TI could develop 
Kelvin-Helmholtz secondary instabilities if cold clouds 
move at transonic speeds relative to a warm background. 
They speculate that such motions could be the result of 
buoyancy forces or some pressure imbalance. However, 
this scenario does not seem to apply to our simulations 
where pressure imbalances become quite small at late 
times. A related possibility would be secondary insta- 
bilities caused by differential rotation. Again, in the 
present simulations this did not occur either. Instead, 
shear mainly causes the flow to become two-dimensional, 
that is, uniform in the streamwise direction. However, in 
the simulations the TI shows no tendency of subsequent 
fragmentation of structures in the streamwise direction. 
There might still be some hope that the fragments could 
be susceptible to a baroclinic instability, but this may re- 
quire substantially higher resolution than what we have 
considered in the present paper. 

In the pure TI, cases the system develops into a new 
segregated state in which each phase is stable. The cold 
patches have a tendency to coalesce into bigger ones that 
are more resistant to the possibility of breaking up. It 
is conceivable that the process of coalescence is slowed 
down when the value of x is decreased. This might be- 
come more plausible when realizing that, because of the 
thermal instability, the energy equation is essentially of 
the type of a reaction-diffusion equation. Under the as- 
sumption of perfect pressure equilibrium at all times, El- 
phick et al. (1991) showed that this equation permits 
traveling kink solutions. If a front were to travel into 
an unstable equilibrium state, the front speed would be 
proportional to the square root of the product of dif- 
fusivity and the growth rate of the instability. In the 
present case, however, warm and cold equilibrium states 
"compete" against each other, so fronts would not propa- 
gate. Only in two and three dimensions, where fronts are 
in general curved, they tend to be driven diffusively to- 
ward the direction of the center of curvature; see Shaviv 
& Regev (1994) as well as Brandenburg & Multamaki 
(2004) for similar results in a different context. How- 
ever, the assumption of perfect pressure equilibrium is 
problematic, because then the density is assumed to be 
inversely proportional to the temperature, so mass con- 
servation is generally not obeyed. In our cases there is 
no perfect pressure equilibrium and one may argue that 
the coalescence is primarily the result of individual dense 
spots moving with the flow toward local and global pres- 
sure minima. 

It is in principle possible that the amount of viscous 
heating might suffice to heat the cold patches enough to 
make them unstable again. However, the amount of vis- 
cous heating is insufficient in all the cases that we have 
investigated. Only when an external forcing function is 
added to give the flow a rms velocity of 10-30 km/s does 
the total amount of heating become comparable to F, 
that is, the level of the imposed uniform heating. Obvi- 
ously, we cannot exclude the possibility of Tl-drivcn tur- 
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bulcncc for smaller viscosity, smaller thermal diffusivity, 
or both. It may therefore be useful to revisit this issue 
in future when simulations at higher resolution become 
more affordable. 

Detailed models of the heating and cooling properties 
of the warm and cold components of the ISM have been 
used to calculate the equilibrium curves, which in prac- 
tice predict the range of stable versus unstable densities, 
temperatures and pressures in the ISM (see, e.g., Wolfire 
et al. 1995, 2003). Calculations, such as those presented 
in this paper, of the onset and nonlinear stages of the 
TI, are needed in order to investigate the actual equilib- 
rium pressure realized in a system described by a certain 
equilibrium curve; from this, the characteristic densities, 
temperatures and filling factors for the warm and cold 
phases can also be determined. 

Our one-dimensional calculations (Figure 5) of the 
standard model of Wolfire et al. (1995) show that the 
mean pressure realized in the unstable regime remains 
roughly constant at 1750 K cm""^ over the whole range of 
unstable densities, {p) = (1-5) x 10~^^gcm~'^, and that 
the pressure is close to the minimal value of 1540 K cm~^. 
We note, however, that the equilibrium curve differs from 
the original one because we have used p,—Q.(!)2 instead of 
/i=l, and it also differs somewhat from Sanchez-Salcedo 
et al. (2002) as we have used revised coefhcients based 



on more accurate continuity considerations. The cor- 
responding temperatures and number densities at this 
equilibrium pressure are Tq = 126 K, nc = 8.6 cm~^ and 
Tw = 9430 K, nw = 0.1 cm^"^. This behavior carries 
over into the three-dimensional regime. The calculations 
with forced turbulence, where the strongest forcing re- 
sults in turbulent pressures exceeding the thermal pres- 
sure by a factor of < 3 show that the mean pressure 
increases only by about 25%, even though the level of 
turbulence is relatively strong. The mean pressure ob- 
tained in this case in the three-dimensional calculations 
is roughly 2170 Kcm^'^, which is in agreement with the 
observed median pressure of p/k^ ~ 2250 K cm""^ from 
Jenkins & Tripp (2001). 
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